Critical Dynamics in a Binary Fluid: Simulations and Finite-size Scaling 
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We report comprehensive simulations of the critical dynamics of a symmetric binary Lennard- 
Jones mixture near its consolute point. The self-diffusion coefficient exhibits no detectable anomaly. 
The data for the shear viscosity and the mutual-diffusion coefficient are fully consistent with the 
asymptotic power laws and amplitudes predicted by renormalization-group and mode-coupling the- 
ories provided finite-size effects and the background contribution to the relevant Onsager coefficient 
are suitably accounted for. This resolves a controversy raised by recent molecular simulations. 
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Introduction. — Thermodynamic and transport 
properties exhibit critical-point singularities with ex- 
ponent values and amplitude ratios that are the same 
for systems belonging to the same universality class. 
As regards static critical behavior, it has been well 
established that fluids of molecules with short-range 
interactions belong to the universality class of three- 
dimensional Ising-type systems [l|. It is expected 
that the dynamic critical behavior of fluids conforms 
to that of model H in the nomenclature of Hohen- 
berg and Halperin 0. Recently, there has been a 
revival of interest in critical phenomena, one reason 
being that computer-simulation techniques have ma- 
tured sufficiently that they can provide interesting de- 
tailed information concerning the static critical behav- 
ior [3|,y,|5|. For instance, recent Monte Carlo simula- 
tions have provided new insights concerning the nature 
of the scaling fields in asymmetric fluids i5| . 

The status of computer simulations of the dynamic 
critical behavior of fluids is much less satisfactory. 
Specifically, on the basis of a recent molecular dy- 
namics simulation of a binary fluid Jagannathan and 
Yethiraj (JY) |^ concluded that the dynamic exponent 
XD that governs the slowing down of critical fluctua- 
tions differs substantially from the value predicted by 
renormalization-group or mode-coupling theory [2, [TJ . 
Sengers and Moldover have pointed out that the con- 
clusion of JY is also in disagreement with reliable ex- 
perimental evidence [8j. 

To address this issue we have undertaken a compre- 
hensive study of the dynamic critical behavior of a sym- 
metric Lennard- Jones mixture (A-l-B) near its conso- 
lute point. We find that the data for the transport 
property that determines the nature of critical slowing- 
down are significantly affected by finite-size effects and 
by a 'background' contribution arising from fluctua- 
tions at small length scales. After properly accounting 
for both these effects our extensive simulations of the 
critical dynamics are fully consistent, both with cur- 



rent theoretical predictions and with the best available 
experimental evidence. 

The model. — Starting from the standard (12, 6) 
Lennard- Jones potential with parameters Eap and 
Cq/? {a, /3 = A, B) we construct a truncated poten- 
tial which is strictly zero for r> re and makes 
both the potential and the force continuous at 
r — Tc Isl- For the parameters, we take 0-^1^1 

= CTBB = CAB = cr, eAA = £BB = '2eAB^£, ^c = 2.5cr, 

and define the reduced temperature as T* = k^T/e. 
The total particle number N = Na+Nq and the vol- 
ume V = L^ are chosen so that the reduced density 
p* = pa^ = N/V is unity and the simulation box edge 
is L/(7 = L* — N^/^. For these parameters the system 
is far from solid-liquid and liquid-gas transitions in the 
temperature regime of interest here. 

In order to evaluate the results of computer simula- 
tions of dynamic critical behavior, accurate information 
regarding the static critical behavior is a prerequisite. 
We have obtained this by using a sen ii-g randcanonical 
Monte Carlo (SGMC) approach H E [ij- In the 
SGMC method, in addition to displacement moves, the 
particles may switch identity (A— >B^A) with both the 
energy change AE and the chemical potential difference 
AiJ, = jjLA—pB entering the Boltzmann factor. For the 
case A/x = 0, of interest here, one has (xa) = (a^e) =1/2 
(with Xa=Na/N) for T > T^. 

Static Properties. — Since our focus is on the dy- 
namic critical behavior, we simply state the results 
found for the static behavior [9|. An accurate, unbiased 
estimate for the reduced critical temperature was ob- 
tained by plotting the fourth-order cummulant Ul — 
((xa-1/2)'')l/((xa-1/2)2)| as a function of T for var- 
ious box sizes L and identifying T(jfrom the asymptot- 
ically common intersection point [a, Isl |l3| : this yields 
T* = 1.4230 ± 0.0005 [i|. The concentration suscep- 
tibility x{T) was calculated from fcsTx = X*T* = 
N{{x%) - (xa)^) (T > Te). The correlation length ^(T) 
was determined from Ornstein-Zernike plots of the 



Fourier transform of the concentration-concentration 
correlation function, Scdq) = T*x*/[1 + q'^S,^ + ...]. In 
the thermodynamic hmit (i.e., in the absence of finite- 
size effects) these properties diverge as x* ~ Toe"''' and 
C ~ £,Q^~'^ where e = {T — Tc)/Tc and we may adopt 
7 — 1.239 and v = 0.629 as the universal critical expo- 
ncnts for the 3-dimcnsional Ising universality class 13^ . 
Our SGMC simulations ^2,) then yield Tq = 0.076±0.006 
and ^o/o- = 0.395 ±0.025. 

Dynamics. — We investigated the dynamic behavior 
by implementing a microcanonical Molecular Dynamics 
(MD) simulation [lj|. For this study, multiple indepen- 
dent initial configurations were prepared from SGMC 
runs with 5x10^ Monte Carlo steps (MCS) per parti- 
cle. This was followed by a microcanonical thermal- 
ization for 2x10^ MD steps in the NVT ensemble us- 
ing the Andersen thermostat [lj| before the production 
runs commenced. For the MD simulations, the particle 
masses were taken equal: niA = "niB = 'ti. The stan- 
dard Verlet velocity algorithm [lj| was employed with 
a time step Jt* =0.01/\/48 [in units Iq = [ma'^ /ef/'^]. 

Self-difFusion. — Restricting attention to temper- 
atures T >Tc and to the critical concentration 
Xc = XA =2:b = 1/2, the symmetry of our model dic- 
tates that the self-diffusion constant is the same for A 
and B particles: Da = Db = D. We have calculated the 
reduced self-diffusion coefficient D* from mean square 
displacements via {a'^/to)D* = D = hmi^oo([?^i(0) — 
ri{t)]'^)/6t, where the average (•) includes all A and B 
particles. The results are shown in Fig. ^a.) as a func- 
tion of e. No anomalous critical behavior is detected 
and the linear behavior is consistent with previous sim- 
ulation studies [y, ll5|. An MD calculation [ig has 
suggested a weak singularity in the self-diffusion near 
vapor-liquid criticality but no corresponding anomaly 
has yet been detected experimentally |l7^ . 

Shear viscosity. — The shear viscosity is expected 
to diverge as S,^" with x„ ^ 0.0679 according to recent 
theoretical calculations [l3| ; this value is in good agree- 
ment with the best available experimental information 
[^ J^j • We have calculated the reduced shear viscosity 
ry* from the appropriate Green-Kubo formula 20] 



77* = {4/aVm'T*) / dt{a,yiO)a,y{t)), (1) 



where the pressure tensor is (Jxy{t) ~ 'Ylii=iV'^i'^ixViy + 
\ ^ji^i) l^« ~ ^jWv{\'''i ~ ^il)] with F and vi the force 
between particles i and j and the velocity of particle 
i, respectively. The numerical data for if obtained 
from simulations with N = 6400 particles are shown 
in Fig. njb). As always in MD simulations, accurate 
estimation of the shear viscosity is difficult and the 
±5% error bars prevent us making any strong state- 
ments about the singular behavior of 77* . But the slight 



increase of 77 as T ^ Tc is actually consistent with 
the predicted power-law divergence 77* « 7706"''^'' with 
v = 0.629 and Xr, = 0.068. The corresponding least- 
squares fit in Fig.^b) yields 770 = 3.87 ± 0.3. 
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FIG. 1: (a) Reduced self-diffusion constant D* for a system 
oi N — 6400 particles as a function of T. The dashed line 
guides the eye and shows that D{Tc) is nonzero, (b) Log-log 
plot of the reduced shear viscosity for the same system. 







I I I I I I I I I 



T 1 1 1 1 1 1 r 



\^oc(T) 



N = 6400 



7K^ I 



1 






I I I I I I I I I 



e=(T-T^)/T^ 

J I I I L 



0.05 0.1 0.3 0.5 



0.8 



FIG. 2: Variation with temperature of the Onsager coef- 
fecient C{T) for a system of A'^ = 6400 particles. Note the 
expansion of the scale for e < 0.1. The dotted line repre- 
sents an effective background contribution: see text. 

Mutual diffusion. — Dynamic renormalization-gro- 
up and mode-coupling theories predict that the mutual 
diffusion coefficients Dab (T) will vanish asymptotically 
as ^~^" , where 



xd 



1 



1.068, 



(2) 



so that there is only one independent exponent charac- 
terizing the dynamic critical behavior of fluids . This 
relation has been verified experimentally |2l|. 

The mutual diffusion coefficient -Dab = ('''^/^o)-Dab 
is related to a corresponding Onsager coefficient £ via 
D^B = C/x* 22] . We have calculated D^^ by adopt- 
ing the result x*(T') « Foe"^ previously obtained, and 



using MD simulations to determine C{T) from the ap- 
propriate Green-Kubo formula [20j 

/•OO 

C{T) = {to/NT*a') dt{J^^{0)J^^{t)), (3) 



^A^B 



where J^^it) = a;B Ef^i ^.,A(i) " ^A EfJl ^»,b W, m 
which Vi^a is the velocity of particle i of species a. 

If, somewhat naively, one fits the numerical values for 
Dab obtained for N — 6400 particles and e > 0.01 to a 
power law of the form i^^B ^ ^^^"^'f one finds a value of 
about 1.6 for the effective critical exponent; this is even 
larger than the corresponding value Xr^ = 1.26 ± 0.08 
derived by Jagannathan and Yethiraj [g from their MD 
simulations! Both values differ substantially from the 
theoretical prediction recorded in |(5J. 

To resolve this issue we must focus on the Onsager 
coefficient C since the simulation data for _y* in our 
model are consistent with Ising criticality JS|. While 
the divergence of x* near Tc is strongly dominated by 
long-range fluctuations, it is known that the Onsager 
coefficient of fluid mixtures near a consolute point (or, 
its equivalent, the thermal conductivity of a fluid near a 
vapor-liquid critical point) contains a critical enhance- 
ment AC{T) due to long-range fluctuations together 
with a significant background which arises from fluc- 
tuations at small length scales [22, 123 and has weak 
temperature dependence l24l |: thus we write 



CiT)^ACiT) + CbiT). 



(4) 



Such a separation has proved essential in reconciling 
experimental data for Dab (T) with theory |2J| . 

In Fig. |2] we show a plot vs. e of the numer- 
ical data obtained for C from the simulations with 
N — 6400 particles. The data do indeed suggest the 
presence of a significant background. Theory predicts 
that A£ satisfies a Stokes- Einstein relation of the form 
AC = i?£)r*x*cr/67r?7*^, where Rd is a universal dy- 
namic amplitude ratio that is of order unity |3, l23l | . It 
thus follows that AC should diverge as 



AC w QT*e- 



with ux — x\v ~ 0.567, (5) 



while x\ = {2L/v) — xd — 0.902. Adopting the value 
Rd ~ 1.05 |23| we find, using the values for Fq, ^o and 
rjo reported above, that a sound theoretical estimate of 
the amplitude Q for our model is 



Q= (2.8 ±0.4) X 10" 



(6) 



Finite-size scaling. — Since the background Cb{T) 
derives from atomic length scales, it should vary little 
with L. However, the possibility of significant finite- 
size effects on the critical part, AC{T), must be rec- 
ognized and allowed for. Note, in particular, that al- 
though static properties may (as here |9j) exhibit negli- 
gible finite-size deviations for the range of (T — Tc) and 
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FIG. 3: Finite-size scaling plots of the critical part of the 
Onsager coefficient for trial values of the effective back- 
ground L\ with y — L/^. We have accepted xx — 0.90 
and set j/o = 7. The filled symbols represent data for dif- 
ferent system sizes at T* = 1.48. The arrow on the right 
marks the theoretical value ||SJ for the critical amplitude Q. 
The dotted lines guide the eye; the dashed line is a scaling- 
function fit embodying the optimal value of C^ . see text. 



L simulated (see Fig.O, the same need not be true for 
transport coefficients. To tackle this problem we write 
the finite-size scaling ansatz 5j 25^ ^6J as 



AC/T* « QWiy)e- 



(7) 



where y = L/^ while W{y) is a finite-size scaling func- 
tion that must vary as Woy^^ [1 -I- Q( ji/^/' ^ )1 f or small y, 
since A£(Tc; L) is finite for i < oo [llillil. For large 
y one may quite generally write 



W{y) = 1 + Wooe-'^yy^ + 



(8) 



where IF(oo)=l is needed to reproduce © when L— >oo 
while, for static properties in short-range systems, n is 
a small integer |5l l25l l26|. However, for dynamic co- 
efficients, where long-time tails, etc., may enter, one 
must be prepared for n = Q implying only an L^"^ de- 
cay of finite-size deviations; the exponent ip demands 
more detailed, currently unavailable theory. 

To analyze the C{T; L) data a scaling plot of 
>Vl(T) = (A£/r*)e'^^ vs. y is desirable: by Q and © 
this should approach Q for large y. But the background 
Cb{T), albeit slowly varying, is unknown! To meet this 
challenge, we introduce an effective background param- 
eter C^ , and adjust it to optimize data collapse: See 
Fig. O which presents Wl (T) for three illustrative val- 
ues of Cl vs. the bounded variable [y/{yo + y)]^^ 



in which, purely for convenience, we have set yo = 7. 
The optimal value, which serves as a rough estimate 
of CbiTc), proves to be Cl^^ ^{3.3 ± 0.8)xlO-3 Q. 
For this assignment we find that a good fit (see the 
dashed hne in Fig. |3Jl is provided by Wl — Q/[l + 
Po/yil + v'^/pDY^ with po = 5.8 ± 0.5 and pi~13.8 
while Q = (2.7±0.4) X 10-3. This estimate for Q is 
in gratifying agreement with the theoretical value re- 
ported in ©. 

The quantitative significance of the finite-size ef- 
fects can be appreciated from Fig. ^ where the scaling- 
function fit has been used to estimate C{T) for N = 
2.56x10^ and N^oo. Note also that the fit for Whiy) 
corresponds to n = and ■0 = 3 in ((S)). Further ex- 
ploration suggests that if an ultimate exponential de- 
cay does arise [if n = f in ©] it sets in only for 
y = L/^ » 30. 

In summary. — The extensive simulations we have 
performed for the transport properties near the demix- 
ing point of our symmetric but otherwise not unreal- 
istic binary fluid model are, when appropriately an- 
alyzed with due attention to strong finite-size effects 
and a background contribution, completely consistent 
with current theoretical predictions and the best avail- 
able experimental data. Not only are exponent values 
and the dynamic exponent relation ((21 respected but 
the amplitude value 10 is also confirmed. While in- 
creased computer power and more refined data analysis 
might eventually provide more stringent tests of theory, 
such as the value of Rd , the necessary resources appear 
rather demanding. 
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FIG. 4: Variation of the Onsager coefficient with T for sys- 
tems of increasing size based on the scaling function fit and 
optimal value of CI [~ Ch{Tc)]: see the dotted line. For 
T within 1% of Tc, reliable estimates of jC{T) would require 
N >3x 10'' and system sizes exceeding 30a. 
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